/*
 * com.javable.cese -- Cell Electrophysiology Simulation Environment (CESE)
 *
 * Copyright (C) 1999 - 2005 Sergey Missan
 * 
 * Dalhousie Univ., Halifax, NS, Canada
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public
 * License along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 * Different license can be applied to this file 
 * under specific agreement with the code author
 * 
 */

package com.javable.cese.templates;

import java.io.Serializable;

/**
 * This class is an abstract template for all AP models.
 * <p>
 * It provides the following basic functionality:
 * <ul>
 * <li>calculates experiment time and stimulation parameters,
 * <li>carries out the main simulation loop,
 * <li>and connects to the <code>AbstractDataExport</code> for exporting data
 * from the model to a persistent storage.
 * </ul>
 * This class is designed to be threadsafe.
 * 
 * @see AbstractDataExport
 */
public abstract class AbstractAPModel extends Object implements Serializable,
        Runnable {

    /**
     * This is a special class for exporting data from the model to a storage
     * class. You should use your own implementation of the
     * <code>AbstractDataExport</code>.
     */
    protected transient AbstractDataExport exporter;

    /**
     * Protection flag for simulation time
     */
    protected boolean isSimulation = false;

    /* Time Step */

    /**
     * Total Experiment Time (ms)
     */
    protected double etime = 800;

    /**
     * Sets Total Experiment Time (ms)
     * 
     * @param setime new value
     */
    public synchronized void setetime(double setime) {
        if (!isSimulation) etime = setime;
        setSteps();
    }

    /**
     * Returns Total Experiment Time (ms)
     * 
     * @return value
     */
    public double getetime() {
        return etime;
    }

    /**
     * Time step (ms)
     */
    protected double dt = 0.05;

    /**
     * Sets Time step (ms)
     * 
     * @param sdt new value
     */
    public synchronized void setdt(double sdt) {
        if (!isSimulation) dt = sdt;
        setSteps();
    }

    /**
     * Returns Time step (ms)
     * 
     * @return value
     */
    public double getdt() {
        return dt;
    }

    /**
     * Initial Time (ms)
     */
    protected double t = 0;

    /**
     * Sets Initial Time (ms)
     * 
     * @param st new value
     */
    public synchronized void sett(double st) {
        if (!isSimulation) t = st;
    }

    /**
     * Returns Initial Time (ms)
     * 
     * @return value
     */
    public double gett() {
        return t;
    }

    /**
     * Number of Steps
     */
    protected int steps = 0;

    /**
     * Sets Number of Steps
     * 
     * @param ssteps new value
     */
    public synchronized void setsteps(int ssteps) {
        if (!isSimulation) steps = ssteps;
    }

    /**
     * Returns Number of Steps
     * 
     * @return value
     */
    public int getsteps() {
        return steps;
    }

    /**
     * Main time loop control variable
     */
    private int increment;

    /* Voltage */

    /**
     * Membrane voltage (mV)
     */
    protected double v = -86.5;

    /**
     * Sets Membrane voltage (mV)
     * 
     * @param sv new value
     */
    public synchronized void setv(double sv) {
        v = sv;
    }

    /**
     * Returns Membrane voltage (mV)
     * 
     * @return value
     */
    public double getv() {
        return v;
    }

    /**
     * Change in Voltage / Change in Time (mV/ms)
     */
    protected double dvdt;

    /**
     * Returns Change in Voltage / Change in Time (mV/ms)
     * 
     * @return value
     */
    public double getdvdt() {
        return dvdt;
    }

    /**
     * New Voltage (mV)
     */
    protected double vnew;

    /**
     * New dv/dt (mV/ms)
     */
    protected double dvdtnew;

    /**
     * Condition to test for dvdtmax
     */
    protected boolean dvdtmax = true;

    /* Total Current and Stimulus */

    /**
     * Time period during which stimulus is applied (ms)
     */
    protected double stimtime = 10;

    /**
     * Stimulus Amplitude (uA/cm^2)
     */
    protected double st = -80;

    /**
     * Sets Stimulus Amplitude (uA/cm^2)
     * 
     * @param sst new value
     */
    public synchronized void setst(double sst) {
        st = sst;
    }

    /**
     * Returns Stimulus Amplitude (uA/cm^2)
     * 
     * @return value
     */
    public double getst() {
        return st;
    }

    /**
     * Time when Stimulus is Applied (ms)
     */
    protected double tstim = 10;

    /**
     * Sets Time when Stimulus is Applied (ms)
     * 
     * @param ststim new value
     */
    public synchronized void settstim(double ststim) {
        tstim = ststim;
    }

    /**
     * Returns Time when Stimulus is Applied (ms)
     * 
     * @return value
     */
    public double gettstim() {
        return tstim;
    }

    /**
     * Stimulus Duration (ms)
     */
    protected double stl = 0.5;

    /**
     * Sets Stimulus Duration (ms)
     * 
     * @param sstl new value
     */
    public synchronized void setstl(double sstl) {
        stl = sstl;
    }

    /**
     * Returns Stimulus Duration (ms)
     * 
     * @return value
     */
    public double getstl() {
        return stl;
    }

    /**
     * Basic Cycle Length (ms)
     */
    protected double bcl = 400;

    /**
     * Sets Basic Cycle Length (ms)
     * 
     * @param sbcl new value
     */
    public synchronized void setbcl(double sbcl) {
        bcl = sbcl;
    }

    /**
     * Returns Basic Cycle Length (ms)
     * 
     * @return value
     */
    public double getbcl() {
        return bcl;
    }

    /**
     * Total current (uA/cm^2)
     */
    protected double it;

    /**
     * Sets Total current (uA/cm^2)
     * 
     * @param sit new value
     */
    public synchronized void setit(double sit) {
        it = sit;
    }

    /**
     * Returns Total current (uA/cm^2)
     * 
     * @return value
     */
    public double getit() {
        return it;
    }

    /**
     * Default constructor
     */
    public AbstractAPModel() {
        setSteps();
    }

    /**
     * This method returns current <code>AbstractDataExporter</code> for this
     * model.
     * 
     * @return current exporter
     */
    public AbstractDataExport getExporter() {
        return exporter;
    }

    /**
     * This method sets given <code>AbstractDataExporter</code>.<br>
     * In fact, you should provide your own class that extends
     * <code>AbstractDataExporter</code>.
     * 
     * @param ex new exporter
     */
    public void setExporter(AbstractDataExport ex) {
        this.exporter = ex;
    }

    /**
     * Implementation of <code>Runnable</code> interface. <br>
     * We need it for threading. It invokes <code>simulate</code>.
     */
    public void run() {
        simulate();
    }

    /**
     * This method contains the main simulation time loop. <br>
     * It also invokes methods in <code>AbstractDataExport</code>:<br>
     * <code>prepareExport()</code> (before the time loop), <br>
     * <code>doExport()</code> (in every time loop iteration), <br>
     * and <code>finishExport()</code> (after the time loop).
     * 
     * @see AbstractDataExport
     * @throws ArithmeticException when there is a problem with our simulation.
     */
    public synchronized void simulate() throws ArithmeticException {
        try {
            isSimulation = true;
            if (exporter != null) {
                exporter.setModel(this);

                exporter.prepareExport();

                /*
                 * Beginning of Time Loop (this is where we begin computing
                 * currents)
                 */
                for (increment = 0; increment < steps; increment++) {
                    if (Thread.currentThread().isInterrupted()) {
                        exporter.finishExport();
                        isSimulation = false;
                        return;
                    }

                    doSimulate();
                    exporter.doExport(increment);
                }

                exporter.finishExport();
            }
            isSimulation = false;
        } catch (Exception e) {
            isSimulation = false;
            e.printStackTrace();
            throw new ArithmeticException("Simulation error");
        }
    }

    /**
     * This method invokes in every time loop iteration in the
     * <code>simulate()</code> method. <br>
     * Since it's abstract, it should be overrided in your model and contain
     * some useful calculations.
     */
    public abstract void doSimulate();

    /**
     * Calculates number of iterations to run the simulation
     */
    private void setSteps() {
        steps = (int) Math.round(etime / dt);
    }

    // //////////
    // These methods implement step size control by calculating an RK4 solution
    // and
    // an RK5 solution. This is used to set local (truncation) error to within
    // user
    // tolerance and adjust stepsize accordingly. RK4 and RK5 are calculated
    // using
    // Fehlberg's methods.
    // This code is based on Jordan Bauman's implementation.
    // //////////

    private double[] states;

    private double lambda = 100;
    private double tolerance = 1e-3;

    /**
     * Returns truncation tolerance
     * 
     * @return tolerance.
     */
    public double getTolerance() {
        return this.tolerance;
    }

    /**
     * Sets truncation tolerance
     * 
     * @param tolerance new value for tolerance
     */
    public void setTolerance(double tolerance) {
        this.tolerance = tolerance;
    }

    private double sigma = .9;
    private double vvals[][] = new double[100][6];

    // initializing the Fehlberg constants:
    private double[] alpha = { 0, 1 / 4., 3 / 8., 12 / 13., 1., 1 / 2.};
    private double[] gamma = { 25 / 216., 0, 1408 / 2565., 2197 / 4104.,
            -1 / 5., 0};
    private double[] delta = { 16 / 135., 0., 6656 / 12825., 28561 / 56430.,
            -9 / 50., 2 / 55.};
    private double[][] beta = { { 0, 0, 0, 0, 0}, { 1 / 4., 0, 0, 0, 0},
            { 3 / 32., 9 / 32., 0., 0., 0},
            { 1932 / 2197., -7200 / 2197., 7296 / 2197., 0., 0},
            { 439 / 216., -8., 3680 / 513., -845 / 4104., 0},
            { -8 / 27., 2., -3544 / 2565., 1859 / 4104., -11 / 40.}};

    /**
     * Returns states for Markovian process
     * 
     * @return current states
     */
    public double[] getStates() {
        return states;
    }

    /**
     * Sets new states array
     * 
     * @param x new states
     */
    public void setStates(double[] x) {
        states = x;
    }

    /**
     * Main method for Markovian calculations. <br>
     * Should be invoked at every time step.
     */
    public void doMarkovian() {
        // double[] y4 = new double[getStates().length];
        double[] y5 = new double[getStates().length];
        // calculate the v values from t and yn
        calculateVValues(t, getStates());
        // Calculate the RK4 and RK5
        // y4 = integrateRungeKuttaFehlberg(getStates(), gamma);
        y5 = integrateRungeKuttaFehlberg(getStates(), delta);
        // local extrapolation
        setStates(y5);
    }

    /**
     * This routine contains the derivative and needs to be changed for
     * different problems
     * 
     * @param x independent variable (usually time)
     * @param y array of arguments
     * @return resulting array
     */
    public double[] derivative(double x, double[] y) {
        return null;
    }

    /**
     * Calculates local truncation error
     * 
     * @param y4 result of 4th order integration
     * @param y5 result of 5th order integration
     * @return truncation error
     */
    public double calculateDelta(double[] y4, double[] y5) {
        double max = 0;
        double maxcheck = 0;
        // we need to find the ||y^-y||(infinity)
        for (int count = 0; count < y4.length; count++) {
            maxcheck = Math.abs(y5[count] - y4[count]);
            max = Math.max(maxcheck, max);
        }

        return max;
    }

    /**
     * Adjusts stepsize based on a local truncation error
     * 
     * @param delta current truncation error
     * @return new step size
     */
    public double updateStepSize(double delta) {
        double steptemp = 0;
        double stepnext = 0;
        if (delta != 0) {
            steptemp = sigma * dt * Math.pow(tolerance / delta, 1 / 5);
            stepnext = Math.max(steptemp, dt / lambda);
            stepnext = Math.min(steptemp, dt * lambda);
        }

        return stepnext;
    }

    /**
     * This routine calculates the v values needed for Fehlberg's method
     * 
     * @param x independent variable (usually time)
     * @param y array of arguments
     */
    public void calculateVValues(double x, double[] y) {
        double[] column = new double[y.length];
        int yindex = 0;
        int yinner = 0;
        int vindex = 0;
        double[] yshift = new double[y.length];
        int shiftindex;
        // the notation is v[i][j] means vj calculated using y[i]
        // first we calculate v[0] for all y's
        vindex = 0;
        column = derivative(x, y);
        // the answer needs to get read into a column
        for (yindex = 0; yindex < y.length; yindex++) {
            vvals[yindex][vindex] = column[yindex];
        } // now the loop to calculate v1->5
        for (vindex = 1; vindex <= 5; vindex++) { // resetting all the
            // elements
            // of yshift = to y
            for (yinner = 0; yinner < y.length; yinner++) {
                yshift[yinner] = y[yinner];
            } // the yshift=yn+h*sum(j=0...i-1)(b[i][j]*v[i])
            for (shiftindex = 0; shiftindex <= vindex - 1; shiftindex++) {
                // it must work out the shift for the whole vector
                for (yinner = 0; yinner < y.length; yinner++) {
                    yshift[yinner] += dt * beta[vindex][shiftindex]
                            * vvals[yinner][shiftindex];
                }
            } // find the value of the derivative with these points
            column = derivative(x + alpha[vindex] * dt, yshift);
            // put the answer into the v array
            for (yindex = 0; yindex < y.length; yindex++) {
                vvals[yindex][vindex] = column[yindex];
            }
        }
    }

    /**
     * Integrates using 5-th order Runge-Kutta-Fehlberg method
     * 
     * @param y array of arguments
     * @param constants v values for the given step
     * @return current function value
     */
    public double[] integrateRungeKuttaFehlberg(double[] y, double[] constants) {
        int yindex = 0, sumindex = 0;
        double cvsum = 0;
        double[] answer = new double[y.length];
        // calculate the y values
        for (yindex = 0; yindex < y.length; yindex++) {
            cvsum = 0;
            // sums up the vij row with weights
            for (sumindex = 0; sumindex <= 5; sumindex++) {
                cvsum += dt * constants[sumindex] * vvals[yindex][sumindex];
            }

            answer[yindex] = y[yindex] + cvsum;
        }
        return answer;
    }

    /**
     * @return Returns the userName.
     */
    public String getUserName() {
        return "--";
    }

    /**
     * @return Returns the organisation.
     */
    public String getOrganisation() {
        return "--";
    }

    /**
     * @return Returns the serialNumber.
     */
    public String getSerialNumber() {
        return "--";
    }

}